##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(DESeq2) 

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineCounts.FilterLowExpression-MergeMutiSample.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/wgcnv"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

dat_tpm_all <- dat_tpm[,-2]

##########################################################################################
## 构建分组信息
condition <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 2))    
coldata <- data.frame(row.names = colnames(dat_tpm_all)[-1] , condition)

## counts变为整数
database <- round(dat_tpm_all[,-1])
rownames(database) <- dat_tpm_all$gene_id

##########################################################################################
## 构建dds矩阵：就是利用上面的counts.txt文件（即对象database）和分组信息（即对象coldata）构建。
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)

##########################################################################################
## 如有多个组需要比较，建议不要将其两两分开而是一起分析，通过在results时指定contrast对象，获得两两的比较结果，
## 利用DESeq（）函数标准化dds矩阵；
dds1 <- DESeq(dds)    # 将数据标准化，必要步骤！！！

##########################################################################################
## 用于加权神经共表达网络
## https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html
## https://www.biostars.org/p/280650/
## We then recommend a variance-stabilizing transformation. 
## For example, package DESeq2 implements the function varianceStabilizingTransformation which we have found useful, 
## but one could also start with normalized counts (or RPKM/FPKM data) and log-transform them using log2(x+1). 
## For highly expressed features, the differences between full variance stabilization and a simple log transformation are small.

dat_forwgcna <- varianceStabilizingTransformation(dds1, blind = TRUE, fitType = "parametric")
dat_forwgcna_out <- data.frame(assay(dat_forwgcna))
dat_forwgcna_out$gene_id <- rownames(dat_forwgcna_out)

image_name <- paste0( out_path , "/CombineCounts.FilterLowExpression-MergeMutiSample.varianceStabilizingTransformation.tsv" )
write.table( dat_forwgcna_out , image_name , row.names = F , sep = "\t" , quote = F )
